function LaplaceJacobigeneral
% Name: LaplaceJacobigeneral.m
% Description: Solves Laplace's equation, d2U/dx2+d2U/dy2=0,
% using Jacobi iteration over a rectangular domain.
% Subfunction: fdistance: computes the distance between 2 N-vectors using the infinity norm.
% Boundary conditions: Dirichlet at top (y max) of domain and zero gradient on
% other 3 sides.
%
% Note: indexing follows the convention that i increases in the x direction and j increases
% in the y direction as per Chapter 4 - this is different to matrix rows and columns.
clear; clc; clf % clear variables, console and figures
%
% Define and initialise problem parameters
xlength=200.0; % length of domain in x direction [m]
ylength=100.0; % length of domain in y direction [m]
%
NX=101; % number of grid points in x direction
NY=51; % number of grid points in y direction
dx=xlength/(NX-1); % grid spacing in x
dy=ylength/(NY-1); % grid spacing in y
x=0:dx:xlength;
y=0:dy:ylength;
%
% Create a matrix for u at computational and ghost grid points. The
% interior unknown u values are initially set to zero. Ghost points
% around left, right and lower sides of the domain are needed to implement
% the zero gradient boundary conditions there.
% Calculation indices are i: 2 to NX+1; j: 2 to NY (data at j=NY+1 is held constant)
grid=zeros(NX,NY); % create initial matrix for grid values
%
% Fill in the **fixed** boundary values at the top (j=NY) of the domain
grid(:,NY)=ylength+0.009*(x.*(xlength-x)); % given formula
% Now create ghost locations around bottom, left and right sides
u0=zeros(NX+2,NY+1);  % correct sized array
% embed grid into u0
u0(2:NX+1,2:NY+1)=grid;
u1=u0; % matrix for updated solution
%
% At the left, right and bottom boundaries
% the gradient of u is zero. Using the central difference the ghost value
% will equal the value interior to the boundary value. Hence the boundary
% values will need to be updated after each iteration.
%
% Parameters for general Jacobi 5-point iteration formula
b=dx/dy;
b2=b*b;
c1=2*(1+b2);
tol=5*10^(-1); % tolerance for the iteration
iterations=0;  % iteration counter
%
disp('program running ...')
% Start Jacobi iteration scheme
dist=tol+1;  % make dist > tol
while(dist>tol)
iterations=iterations+1;
% insert ghost values based on zero gradient boundary conditions
u0(1,:)=u0(3,:);  % left ghost value
u0(NX+2,:)=u0(NX,:); % right ghost value
u0(:,1)=u0(:,3);  % lower ghost value
%
% Jacobi iteration
for j=2:NY;
for i=2:NX+1;
u1(i,j)=(b2*u0(i,j-1)+u0(i-1,j)+u0(i+1,j)+b2*u0(i,j+1))/c1;
end
end
%
% calc distance between solutions at consecutive iterations
% at interior grid points
dist=fdistance(u1(2:NX+1,2:NY),u0(2:NX+1,2:NY));
u0(2:NX+1,2:NY)=u1(2:NX+1,2:NY); % update u values
end % of while loop

% Screen Output
disp('Jacobi iterative scheme has converged')
disp('iterations=')
iterations
disp('see contour plot in graphics window')
[X,Y]=meshgrid(x,y); % for contour plot
[C,h]=contour(X,Y,transpose(u0(2:NX+1,2:NY+1)));
clabel(C,h) % lable contour lines
title('numerical solution')
xlabel('x (m)');
ylabel('y (m)');
disp('program finished')
% end of LaplaceJacobigeneral.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d=fdistance(A,B)
% File name: fdistance.sci
% Description: Function to compute the distance between 2 N-vectors using the infinity norm.
[ma,na]=size(A);
[mb,nb]=size(B);
if ((ma ~= mb)|(na ~= nb) )
disp('matrices are not of equal size, stopping program')
exit
end
diff=A-B;
d=max(max(diff)); % diff is a matrix so need 2 max functions
% end of function fdistance